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Abstract 

The Prokof 'ev Svistunov worm algorithm was originally developed 
for models with nearest neighbor interactions that in a high tempera- 
ture expansion are mapped to systems of closed loops. In this work we 
present the surface worm algorithm (SWA) which is a generalization 
of the worm algorithm concept to abelian Gauge-Higgs models on a 
lattice which can be mapped to systems of surfaces and loops (dual 
representation). Using Gauge-Higgs models with gauge groups Z 3 and 
U(l) we compare the SWA to the conventional approach and to a local 
update in the dual representation. For the Z 3 case we also consider 
finite chemical potential where the conventional representation has a 
sign problem which is overcome in the dual representation. For a wide 
range of parameters we find that the SWA clearly outperforms the 
local update. 



1 Introduction 



Monte Carlo simulations are a powerful tool for the analysis of spin systems 
and lattice field theories and Monte Carlo techniques have seen a tremendous 
development over the last decades. An important aspect of this development 
is the choice of the representation of a physical system that is optimal for 
the Monte Carlo simulation. 

A prominent example for the success of a Monte Carlo simulation in 
an alternative representation is the Prokof'ev Svistunov worm algorithm 
PQ. Originally it was proposed for the simulation of spin systems in a loop 
representation. The loop representation (or dual representation) is obtained 
from the usual spin language by a high temperature expansion where the 
new degrees of freedom are link occupation numbers subject to constraints 
at the sites of the lattice, such that admissible configurations correspond to 
loops on the lattice. The worm algorithm not only solves the problem of 
properly taking into account the constraints in the Monte Carlo update but 
turned out to be outperforming many previous simulation approaches in the 
conventional formulation [2]. 

The worm algorithm concept found many interesting applications also 
for quantum field theories on a lattice. In this area a strong motivation for 
dual representations is the study of quantum field theories with a chemical 
potential, where in many cases the standard representation has complex 
action and a direct Monte Carlo simulation is not possible. Lattice field 
theories that were studied with worm-type algorithms comprise scalar field 
theories HI fermion systems in various settings, in particular with four 
fermi terms or in the strong coupling limit [6], as well as effective theories 
for the QCD phase diagram [7j . All these systems have in common that the 
interaction on the lattice is either supported on a single site or on nearest 
neighbors. The resulting dual representation thus consists of loops. 

A genuinely new element appears in the dual representation of gauge 
theories. There the interaction is based on the plaquettes of the lattice and 
the corresponding dual variables (integers assigned to the plaquettes) form 
surfaces. While for the non-abelian case the structure is rather involved 
[H] , abelian gauge theories have a straightforward representation in terms of 
closed surfaces. Nevertheless only a few suggestions and attempts for a dual 
simulation of abelian gauge theories can be found in the literature [HOE] 
and the main obstacle for a worm-type algorithm is to efficiently generate 
the closed surfaces of the dual representation. 

It is interesting to note that the situation is simplified, when matter is 
coupled to abelian gauge fields: The dual variables of matter fields are fluxes 
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based on the links of the lattice that serve as boundaries of the surfaces rep- 
resenting the gauge degrees of freedom. Despite the fact that an additional 
field appears, the dual representation and in particular its Monte Carlo sim- 
ulation become simpler because the algorithm now also may use plaquettes 
bounded by matter flux. A first analysis of a Gauge-Higgs system in the 
dual representation with a local Monte Carlo update was presented in |10| . 

In this article we now present a new Monte Carlo strategy: The surface 
worm algorithm (SWA) which is a generalization of the worm algorithm to 
a system of surfaces and loops, i.e., dual representations of abelian Gauge- 
Higgs models. The SWA uses two main elements: Changing the flux at an 
individual link as well as changing a plaquette occupation number and the 
flux on two of the links of that plaquette. These steps are used to efficiently 
build up filament-like structures where the link and plaquette occupation 
numbers are altered. We verify and test the surface worm algorithm for 
lattice Gauge-Higgs models with gauge groups U(l) and Z3. The latter 
case has a complex action problem in the conventional approach which is 
overcome by the dual representation. We find that in both models the 
surface worm algorithm outperforms local updates. 

2 Two abelian Gauge-Higgs models and their dual 
representation 

We use two different Gauge-Higgs models based on the gauge groups Z3 and 
U(l) to test the surface worm algorithm and explore its properties. This 
section defines the two models in their conventional representations and 
summarizes their dual form in terms of loops of flux and surfaces. For the 
actual derivation of the dual representation we refer to the literature, to [10] 
for the case of the Z3 model and to [11] for U(l). 

2.1 The Z 3 Gauge-Higgs model 

In the conventional form the degrees of freedom of the Z3 Gauge-Higgs model 
are the gauge fields U xll , v = 1, 2, 3, 4 living on the links of a 4-dimensional 
lattice, and the scalar matter fields <j) x located on the sites. Both sets of 
degrees of freedom are in the gauge group Z3 = {1, e* 27r / 3 , e~* 27r / 3 }. The 
lattice we consider has size V4 = Nf x N t and we use periodic boundary 
conditions for both fields. The action S is a sum of the gauge action Sq and 
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the action Sm f° r the matter fields. The gauge action is given by 

s g = -f EEI^ + u iu P ] > (!) 

where = U X)1/ U x+ p jP U* + ^ v U* tP and /3 is the inverse gauge coupling. 

The action for the matter fields is 

Sm = -kJ2[^ v ' A €U x ,u(I> x +u + e-^^U,;^] , (2) 

x,u 

where a chemical potential fi is coupled to the terms in the temporal direc- 
tion and the hopping parameter k is a positive real number. The partition 
sum of the conventional representation is given by Z = Yl{u<i>} e~ Sa ~ Sa 
where the sum is over all possible field configurations. We stress that in the 
conventional form the Z3 Gauge-Higgs model has a complex action problem 
at non-zero chemical potential, i.e., Sm is complex for /i > 0. 

The partition sum can be rewritten exactly [10J into a dual representation 
where the new degrees of freedom are link variables l xv G { — 1,0, +1} and 
plaquette variables p x ,pu £ {—1, 0, +1}. The partition function is a sum over 
all configurations of the link and plaquette variables, 

Z = C W[p,l]Cs[l]C L [p,l] . (3) 
The configurations {p, 1} in ^ come with real and positive weight factors 

wm = (nn^ pl ) (nrN M ) (n*0 ■ ( 4 ) 

^ x v<p ' ^ x i=l ' ^ x ' 

with 

e 2n _ e - K e P _ p -P/2 



The overall constant C in (ph is given by (3B^B^) Vi . The last contribution 
to the weight Q contains the chemical potential. It is a product over factors 



M ixA with 

Mi = - 
3 



(6) 



where I = +1,0, —1. Note that the factors Mj are real and positive also for 
fi > 0. Thus the complex action problem is solved in the dual representation. 
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The configurations {p, 1} are subjects to the constraints 



x v=l \p:v<p p:v>p 
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c s [i) = U T E^-w ' ( 7 ) 

x \i/=l / 

that both contain the triality function T{n) which is defined to be 1 if n is a 
multiple of 3 and vanishes otherwise. The constraint Cs[l] is a product over 
sites x of the lattice and enforces that the total flux Y^A^x-u,u — lx,v] at the 
site x is a multiple of 3. The constraint Cl[p,1] is a product over links of 
the lattice and forces the combined flux from the plaquettes attached to the 
link and the corresponding link variable to be a multiple of 3. 

The admissible configurations of the dual variables p and / have the 
interpretation of surfaces made of non-zero plaquette variables p x ,up- The 
surfaces can either be closed (without boundaries) or they are bounded by 
loops of link variables that compensate the flux at the links that constitute 
the boundary of the surfaces. 

2.2 The U(l) Gauge-Higgs model 

In the U(l) Gauge-Higgs model the degrees of freedom are gauge fields 
U XfV G U(l) at the links of the lattice and a charged scalar Higgs field 
4> x G C, attached to the sites. Again we consider a 4-dimensional lattice 
with V4 = Ng x Nt and periodic boundary conditions for both fields. The 
gauge action Sg has the same form as in Q - only the link variables are 
U(l)-valued now. The action for the matter fields is given by 



S M = E 



k\4> x \ 2 + x\4) x \ A - ^2 [<t>* x u x>u 4> x +p + <i>xU^ v 4> x+i ) 



(8) 



The parameter k denotes 8 + m 2 , where m is the bare mass parameter and 
A is the quartic coupling. The partition sum Z = f D[U]D[(j)]e~ SG ~ SH is 
given as an integral over all field configurations. Again the partition sum can 
be mapped exactly to a dual representation. Here we need two sets of link 
variables, l XjU G Z, Z Xj „ G Mo, and plaquette occupation numbers p x , P u G Z. 
The dual partition function is a sum over all configurations of the /, / and p 
variables, 
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z = E E ft *] w g w c ^ w Ci [p- z i • ( 9 ) 

{1,0 M 

The weight factors are 

w m[U] = U m 1,7 M7 ■ n p (E0^1+l^l + 2 (W+W)]) ) 

= n > ( io ) 

where / p (/3) denotes the modified Bessel functions and the P(n) are the ele- 
mentary integrals P(n) = dx x n+1 e~ Kx2 ~ Xx4 . In a numerical simulation 
the -P(n) are pre-computed and stored for a sufficient number of values n 
so they can be used for determining the Metropolis acceptance probabilities 
efficiently. Only the I and the p variables are subject to constraints given 
by ( 5{n) is here used to denote the Kronecker delta 5 n fl ) 

cl[p,i\ = n n 5 ( e & x ' u p ~ Px-p,u p \ - ^2 ^ px 'p u ~ Px-p,p»} + ) ' 

x v=\ \ p:v<p p:v>p J 

c s [i] = n^fE^-^-^y ( n ) 

x \u=l J 

The constraints have the same form as for the Z3 case, i.e, we have con- 
straints Cs[l] that are based at the sites for the variables I and constraints 
Cl\p, I] that are based on the links and combine p and I variables. The only 
difference is that the triality functions of ([T]) are for the U(l) case replaced 
by Kronecker deltas, implying that all fluxes must vanish exactly and not 
only modulo 3 as in the Z3 case. 



3 Monte Carlo simulation 

In this section we describe the surface worm algorithm (SWA). We also dis- 
cuss a local Metropolis algorithm (LMA) for the dual representation which 
will be used for cross-checking the results from the SWA. Since the steps 
used in the SWA may be viewed as a decomposition of the local update into 
smaller elements we first discuss the local update. 

For the U(l) Gauge-Higgs model in addition to the plaquette variables 
p and the constrained flux variables I we also have the unconstrained link 
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variables I. Due to the absence of a constraint we can update the link vari- 
ables I using conventional Metropolis techniques, which are well documented 
in textbooks (see, e.g., [12]) and thus are not discussed in this paper. The 
update for the constrained variables discussed here is understood in a back- 
ground configuration of the I variables and in the numerical tests presented 
in Section 4 we simply alternate the update of the constrained variables with 
sweeps for the I fluxes. 

3.1 Local algorithm for the dual representation 

The central aspect of a Monte Carlo simulation in the dual representation 
is to generate only admissible configurations, i.e., configurations that obey 
all constraints. The strategy which we adopt for the local update is to 
start from a configuration where all constraints are obeyed - typically the 
configuration where all flux and plaquette variables are set to - and then to 
offer local changes of the dual variables that do not violate the constraints. 

The simplest local change is to increase or decrease a plaquette occupa- 
tion number p XtUp by ±1 and to compensate the violation of the constraint 
on the links of the lattice by changing the link fluxes l X)(T by ±1. The two 
possible changes (one for increasing p XtUp , one for decreasing) are illustrated 
in Fig. [T] The change of p x ,up by ±1 is indicated by the signs + or — , while 
for the flux variables we use a dashed line to indicate a decrease by —1 and 
a full line for an increase by +1. It is easy to see that the pattern of changes 
for the flux variables not only compensates the violation of the link-based 
constraints from changing p X)Up but also leaves intact the site-based con- 
straints at all four corners of the plaquette. We stress that for the case of 
gauge group Z3 addition of ±1 is understood modulo 3, which is the usual 
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Figure 1: Plaquette update: A plaquette occupation number is changed by 
+1 (lhs. plot) or —1 (rhs.) and the fluxes at the links of the plaquette are 
changed simultaneously. We use a full line for an increase by +1 and a 
dashed line for a decrease by —1. The directions 1 < v < p < 4 indicate the 
plane of the plaquette. 
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Figure 2: Cube update: The plaquette occupation numbers of a 3-cube are 
changed according to the two patterns we show. The edges of the 3-cube 
are parallel to the directions 1<u<p<ct<A. 

addition, except for the cases 1 + 1 = — 1 and —1 — 1 = 1. 

A full sweep of these "plaquette updates" consists of visiting all plaque- 
ttes and offering one of the two changes of Fig. [l] with equal probability. The 
offer is accepted with the usual Metropolis probability min(l, W/ oc /Wi oc ) 
where W/ and Wi oc are the local weights of the trial configuration and 
the old configuration. They can easily be evaluated from the weight factors 
discussed in the previous section. 

It is easy to see that the plaquette update alone is ergodic. Nevertheless 
we found it advantageous to augment the plaquette update with a "cube 
update" that involves only changes of plaquette numbers p x ,up- The cube 
update helps to decorrelate the system in parameter regions where link flux 
has a very small Boltzmann weight. The plaquettes on the faces of 3-cubes 
of our 4-D lattice are changed according to one of the two patterns shown in 
Fig. [2] (for Z3 addition is again modulo 3). The two possibilities are offered 
with equal probability and it is easy to check that the link-based constraints 
are not violated, and since no flux variables are involved also the site-based 
constraints remain intact. A full sweep of cube updates consists of visiting 
all 3-cubes, offering one of the two changes and accepting them with the 
Metropolis probability computed from the local weight factors. 

3.2 Surface worm algorithm 

The surface worm algorithm (SWA) is constructed by breaking up the pla- 
quette update discussed in the previous subsection into smaller building 
blocks used to grow filament-like clusters on which the flux and plaquette 
variables are changed. We will first discuss in detail the SWA for the Z3 
Gauge-Higgs model and then address the modifications necessary for U(l). 
As for any worm algorithm, in the SWA the constraints are temporarily 
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Figure 3: Examples of positive (lhs.) and negative segments (rhs.) in the v- 
p-plane (y < p). The plaquette occupation numbers are changed as indicated 
by the signs. The links marked with full (dashed) lines are changed by +1 
(— 1). The empty link shows where the segment is attached to the worm 
and the dotted link is the new position of the link Ly where the constraints 
are violated. 



violated at a link Ly and the two sites at its endpoints. This is done by 
changing the flux at a randomly chosen link by ±1 (addition is again modulo 
3 for the Z3 case). The defect at Ly is then propagated through the lattice 
by offering steps where a plaquette occupation number is changed by ±1 
and two flux variables at two of the links of the plaquette. We refer to these 
structures as "segments" and show some examples in Fig. |3j Attaching 
segments propagates the link Ly where the constraint is violated through 
the lattice until the worm decides to terminate with the insertion of another 
unit of link flux. Each step is accepted with a Metropolis decision. 

Fig. [3] shows some examples of segments that are used by the SWA. The 
plaquette occupation numbers are changed by ±1 as indicated and also the 
fluxes at two of the links of the plaquette (again we use a full line if the flux 
at a link is increased by +1 and dashed lines for a decrease by —1). We refer 
to a segment as a "positive segment" if the plaquette occupation number 
is increased (first and second segment shown in Fig. [5J) and use "negative 
segment" otherwise (third and fourth segment). The empty link represents 
the link where a segment is attached to the existing filament-like structure 
of the SWA and the dotted link is the new (= shifted) position of the link 
Ly where the constraints are violated ( "head of the worm" ) . 



Thus the SWA proceeds as follows (for an example see Figs. 4(a) -4(e)): 



The SWA starts at a randomly chosen link Lq where the flux is changed 
by ±1 (in the example Fig. |4(a) the flux is changed by — f ). At this 



link and at its endpoints the constraints are violated, i.e., Ly = Lq. 
• Subsequently the SWA either moves Ly by attaching a suitable seg- 



ment (Figs. 4(a) -4(d) ) or decides to change the link flux at Ly to heal 



the violated constraint thus terminating the worm (Fig. 4(e)). 



S 




V 

(a) The worm starts by decreasing the flux in v direction. Subsequently it 
adds a segment in the v-p plane. 




(c) The worm adds a segment in the v-p plane. 




(d) A segment in the v-a plane is added. 




(e) The worm decides to saturate the violated constraint by adding a unit 
of flux in the a direction and terminates. 



Figure 4: Example of a surface worm algorithm on an initially empty lattice. 
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Figure 5: This figure depicts the constraints of the dual partition function. 
It can be used to determine whether a positive or negative segment will be 
inserted by the worm. See the text for a more detailed explanation. 



Whenever the worm decides to add a new segment it first randomly deter- 
mines a new plane for the segment. This plane has to contain the direction 
of the link Ly that currently violates the constraint. Subsequently the worm 
has to determine whether to insert a positive or a negative segment to create 
only admissible configurations. The following steps and Fig. [5] explain how 
the worm selects an admissible segment (l<^<p<<7<4): 

1. Depending on the direction of the link Ly (i.e., Ly || u, Ly \\ p or 
Ly || a) identify Ly as the central link surrounded by four plaquettes 
in one of the diagrams of Fig. [5| 

2. Identify the "old plane" and the "new plane": 

Old plane: plane of the last successfully updated segment. 
New plane: plane of the new trial segment. 

3. If the plaquettes in the old and the new plane are marked by different 
lines (full versus dashed) keep the same type of segment. Otherwise 
change the type of segment from positive to negative or vice-versa. 

Note that when the worm attempts to revisit the last updated plaquette 
(i.e., it moves backwards) then the new and old planes coincide. Thus the 
segment changes and the last move of the worm is undone. 

In addition to the example of Fig. |4j in Fig. [6] we show a short worm 
that generates the plaquette update of the local algorithm discussed in the 
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(a) The worm starts by increasing the flux in the ^-direction and then 
adds a segment in the v-p plane. 




v 

(b) The worm decides to saturate the violated constraint by de- 
creasing the flux at Ly by one unit and terminates. 



Figure 6: Example how the worm generates the local plaquette update dis- 
cussed in the previous subsection. 

previous subsection. We have already stressed that the plaquette update is 
ergodic and the SWA thus is ergodic too. 

The pseudo-code listed below describes the algorithm. For the coordi- 
nates of plaquettes we use P, and L for the coordinates of links. In particular 
the link where the constraints are violated (head of the worm) is denoted 
by Ly. By sp = (pp,lo,h) we denote the current occupation numbers of 
a segment, i.e., the occupation number pp of the plaquette at P and the 
two links fluxes lo and l\ which are changed in the type of segment chosen 
(in the examples of segments shown in Fig. [3] Iq and l\ are the link fluxes 
marked by full or dashed lines). The variable A s = (S p , 5i , 5i x ) denotes the 
change of the occupation numbers of sp. Note that the sign of the change 
AS ( "positive segment" versus "negative segment" ) has to be chosen accord- 
ing to the rules stated in the discussion of Fig. [5j By x © y we denote the 
addition modulo 3 which is the usual addition operation except in the cases 
+1 © +1 = —1 and —1 © —1 = +1. By weight_ratio(6 a) we denote the 
ratio yVi' oc /yVi oc when changing an element a into b. Here a and b are either 
a link flux before and after the change by ±1 or a full segment (a plaquette 
number pp and two link fluxes Iq, Zi) before and after the respective changes. 
Finally, rand() is a random number generator for uniformly distributed real 
numbers in the interval [0, 1). 
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Pseudocode for surface worms: 



select a lattice link Lq randomly 
select 5i G {— 1,+1} randomly 

V <— h Si 

if randO < weight _rat io (/' <— II) 

Il <— 

Ly <S Lo 

else 

terminate worm 

end if 

repeat until worm is complete: 
select a direction p G {±1, ±2, ±3, ±4} 
if p || Ly then 

select (5/ such that the violated constraint at Ly is healed 

V < — h v e 5i 

if randO < weight.ratio (/' ^— Il v ) 

h v < — I' 

terminate worm 
end if 

else 

the plaquette P for a new segment is spanned by Ly and p 
select Ly ^ Ly randomly among the links of P 
choose A s such that the constraint at Ly is healed 

s' < — spffi A s 

if randO < weight_ratio (s' ^— sp) 

sp < — s' 

Ly i — Ly 
end if 

end repeat until worm is complete 

It is straightforward to show detailed balance using the Boltzmann weights 
and that the algorithm is ergodic. 
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Modifications for the U(l) gauge-Higgs surface worm algorithm 



From Eq. ^ and Eq. ^ we observe that the SWA has to be adapted in order 
to simulate the U(l) model: 

• Due to the extra unconstrained set of link variables l x ,u, for U(l) a full 
sweep consists of a worm sweep (V4 = N^Nt worms) to update the l x>v 
and p x ,up plus a conventional local Metropolis sweep to update all l x ,v 

• To extend the range of the constrained variables to all integer numbers 
and enforce the total flux at every link and site to vanish, the operation 
x@y is replaced by a normal addition x + y. In the pseudo-code: sp © A s 
is replaced by sp + A s . 



4 Assessment of the surface worm algorithm 
4.1 Validity of the SWA 

To evaluate the validity of the algorithm we will use several thermodynami- 
cal observables and their susceptibilities. For both models we study the first 
and second derivatives with respect to the inverse gauge coupling /3, i.e., the 
plaquette expectation value and its susceptibility, 

{U) = 6WfN t §p lU Z ' XU = WfN t W lnZ - (12) 

For the Z3 case we also consider the particle number density n and its sus- 
ceptibility which are the derivatives with respect to the chemical potential, 

Id Id 2 

n = — 5 — „ In Z , y„ = — 5 — — -k In Z . (13) 

Finally, for the U(l) model we analyze the derivatives with respect to k, 

wv=imi lnZ • ™ = iw,& XnZ - (14) 

The correctness of the flux representation has already been established 
in |10l ITT] . Thus here we can focus on the SWA. To check for correctness 
we compare the SWA results to the data coming from the local Metropolis 
algorithm (LMA) in the flux representation and for the cases where there is 
no sign problem also to results from a conventional approach in the standard 
representation. 
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Figure 7: Z3 model: (U) and \u at k = 0.5 and fi = as a function of (3 on 
a 10 4 lattice. We compare the results of the SWA (asterisks) to the LMA 
(circles) and the conventional approach (crosses). 



For all simulations we used thermalization and decorrelation sweeps (see 
below for their numbers). For the SWA one sweep consists of V4 = N^Nt 
worms and for the case of U(l) also of a sweep through all unconstrained link 
variables l x ,v For the LMA a sweep is defined as a sequence of plaquette 
updates for all 6V4 plaquettes plus cube updates for all AV4 cubes. For 
the U(l) model and the Z3 case at \i = we can also compare to the 
conventional approach where as usual a sweep is defined as applying one 
local Metropolis update to all degrees of freedom. All error bars we show 
were determined using a Jackknife analysis and are corrected with the factors 
from the respective autocorrelation times (see below). 

For the Z3 model we compared simulations for several parameter sets 
and found very good agreement of the results from the different approaches. 
As examples we show results for two parameter sets: 1) The behavior across 
a crossover transition as a function of /3 at k = 0.5 and \i = (no complex 
action problem) on a 10 4 lattice (Fig. [7}. 2) The behavior across a first order 
transition as a function of [i at k = 0.1 and (3 = 0.6 on a 8 3 x 50 lattice 
(Fig. [8]). In the latter case the standard representation has a complex action 
problem and we only can compare the results from SWA and LMA. For both 
tests we used 10 6 equilibration steps and 10 6 measurements separated by 10 
steps for decorrelation. 

Similarly we also confirmed the correctness of the SWA in the U(l) model 
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Figure 8: Z3 model: The observables (U), xu-, n and Xn as a function of /x 
at k = 0.1 and (3 = 0.6 on a 8 3 x 50 lattice. We compare the results from 
the SWA (asterisks) and the LMA (circles). 



checking the agreement of all three approaches at different parameters and 
lattice sizes. As an example, Fig. [9]shows the results obtained with the LMA 
(crosses), with the SWA (circles) and the conventional approach (asterisks) 
at A = 1 and k = 5, 8 and 9 on a 10 4 lattice. For this test we used 
10 5 equilibration steps and 10 5 measurements separated by 10 steps for 
decorrelation. As for the Z3 case we find very good agreement among the 
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Figure 9: U(l) model: Observables as a function of f3 at A = 1.0 for k = 5, 8 
and 9 on a 10 4 lattice. We compare results from three algorithms: The 
conventional approach (asterisks), the SWA (circles) and the LMA (crosses). 



different approaches thus establishing the correctness of the SWA also for 
the U(l) model. 

4.2 Characteristic quantities of the algorithms 

For a meaningful comparison of the performance and autocorrelation times 
of the SWA and LMA algorithms we study suitable characteristic quantities 
in order to describe the behavior of both algorithms in different regions of 
the parameter space. The definitions are patterned after related quantities 
introduced for the analysis of worm algorithms with open ends [7J. 
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• Plaquette changes V: 



V = average number of plaquettes changed per update 

• Starting fraction S: 

number of successful update starts 

S = < 1 

number of all start attempts 

• Cost ratio C: 

number of attempted changes (plaquettes and links) ^ 
number of accepted changes ~~ 

In these definitions "update" refers to one surface worm for the SWA case. 
For the LMA it is the average of a plaquette and a cube update which we 
consider in a mix of 6V4 plaquette updates and 4V4 cube updates per LMA 
sweep (see above). From the definition of these characteristic quantities it 
is obvious that an optimal algorithm is characterized by a large value of V 
and values of S and 7Z close to 1. 

In Table [T] we show the characteristic quantities for the SWA and LMA 
algorithms in the Z3 case. We compare three different sets of parameters 
denoted by Z-l, Z-2 and Z-3 (see the first column for the corresponding 
parameter values) and four different volumes (second column). The param- 
eters of Z-l are located below the condensation transition shown in Fig. [8j 
the set Z-2 is in the condensed phase (compare Fig. 11 from [TU]) and the 
set Z-3 is inside the crossover region of Fig. [7j 

Table [T] demonstrates that the SWA has a larger probability for starting 
an update than the LMA (Sswa > Slma f° r all data sets and volumes). 
Furthermore the cost ratio 1Z of the SWA is smaller or equal (equal only 
for the set Z-2) to the LMA case. These two quantities indicate that the 
SWA is more effective than the LMA. The observation that V is larger for 
the LMA is mainly due to the fact that an accepted cube update of the 
LMA changes 6 plaquettes (although at the cost of a low acceptance rate). 
It is interesting to note that the values for the characteristic quantities are 
essentially independent of the volume. 

Table [2] collects the data for the U(l) case. Here we consider three 
different sets of parameters U-l, U-2, U-3 (first column) on four different 
volumes (second column). The set U-2 is located very close to the transition 
shown in Fig. [9j the set U-l is below and the set U-3 above the transition. 

The general behavior for the characteristic quantities is essentially the 
same as in the Z3 case: For all sets the starting probability of the SWA 
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Parameters 


T/ 


SsWA 


'PsWA 


CsWA 


Slma 


1 LMA 


Clma 


oet: Zj-1 


4 X OU 


U.zUo 


U.UyO 


b.oyz 


z.ye-o 


A A Kd 

4.400 


ozU. / 


K — U.l, 


o X OU 


U.zUo 


U.UyO 


o.oyz 


z.ye-o 


A AKK 
4.400 


ozU. 1 


8 — n fi 


19 3 x 50 


903 


095 


6 899 


9 9p-3 


4 455 


390 7 


// = 2.0 


16 3 x 50 


0.203 


0.095 


6.892 


2.9e-3 


4.455 


320.7 


Set: Z-jg 


4 3 x 50 


0.245 


1.196 


5.319 


0.172 


5.384 


5.346 


re = 0.1, 


8 3 x 50 


0.244 


1.186 


5.431 


0.172 


5.384 


5.346 


/3 = 0.8, 


12 3 x 50 


0.245 


1.199 


5.320 


0.172 


5.384 


5.346 


/i = 1.6 


16 3 x 50 


0.244 


1.187 


5.425 


0.172 


5.384 


5.346 


Set: Z-3 


4 4 


0.697 


0.802 


3.081 


0.098 


1.286 


10.88 


re = 0.5, 


8 4 


0.698 


0.802 


3.081 


0.098 


1.286 


10.88 


P = 0.28, 


12 4 


0.698 


0.802 


3.081 


0.098 


1.286 


10.88 


/i = 0.0 


16 4 


0.697 


0.802 


3.081 


0.098 


1.286 


10.88 



Table 1: Characteristic quantities for the Z3 model (see the text for their 
definitions). We used 10 6 steps for equilibration and 10 6 measurements 
separated by 2 steps for decorrelation. The errors are smaller than the last 
digit we show. 



is larger than that of the LMA, and also the cost efficiency is considerably 
better for the SWA. As in the Z3 case we find that the average number 
of updated plaquettes V is larger for the LMA, which also here is due to 
the cube updates, which, however, have a much lower acceptance rate as is 
obvious from S and C. The difference in the characteristic quantities between 
the 4 and larger volumes for the set U-2 is due to finite-size effects: In 
the smallest volume the transition is rounded and slightly shifted towards 
smaller values of /3, such that for the smallest volume the parameters we 
work at are further remote from the transition and both algorithms are 
more efficient. 

Finally, comparing Rswa and Sswa for both the Z3 and U(l) cases, 
we observe that even though many worms start successfully, not all of them 
create non-trivial changes, i.e., there is a sizable probability that in the 



second step a worm reverts its initial step. This is also reflected in Fig. 10 



where we show the abundance distribution of the worms as a function of 
their length I defined as the number of segments of a worm. The distribution 
decreases roughly exponentially with I. However, as we shall see in the next 
subsection, a few long worms are enough to have a very efficient sampling. 
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Parameters 


T/ 


SsWA 


TsWA 


CsWA 


Slma 


1 LMA 


Clma 


Set : U-1 


4 


n om 
U.ZU1 


n not; 


o.oyy 


l.ze-o 


1 077 

l.z/ / 


Clf\A P. 


K = 5, 


a4 
o 


n om 
U.ZU1 


n hqc: 
U.Uoo 


a nno 


l.ze-o 


l.z r o 


nno i 

yuy.z 


A — 1 

A — i, 


in4 

XZi 


D 9D1 


U .uou 


6 QD9 


1 9p-3 


1 978 


qnq 4 

i/Uf .4: 


P = 0.40 


16 4 


0.201 


0.085 


6.902 


1.2e-3 


1.278 


909.4 


Set: 17-jg 


4 4 


0.681 


1.275 


3.310 


0.167 


1.813 


6.263 


K = 5, 


8 4 


0.220 


0.199 


6.124 


4.6e-3 


2.243 


224.3 


A = l, 


12 4 


0.220 


0.198 


6.124 


4.6e-3 


2.243 


224.3 


/3 = 0.65 


16 4 


0.220 


0.198 


6.124 


4.6e-3 


2.243 


224.3 


Set: 17-3 


4 4 


0.107 


0.100 


8.775 


0.061 


5.962 


14.82 


« = 8, 


8 4 


0.107 


0.100 


8.773 


0.061 


5.962 


14.92 


A = l, 


12 4 


0.107 


0.100 


8.774 


0.060 


5.962 


14.91 


/? = 1.10 


16 4 


0.107 


0.101 


8.766 


0.060 


5.962 


14.91 



Table 2: Characteristic quantities for the U(l) model (see the text for their 
definitions). We used 10 6 steps for equilibration and 10 6 measurements 
separated by 2 steps for decorrelation. The errors are smaller than the last 
digit we show. 

4.3 Autocorrelation times 

In this subsection we analyze the integrated autocorrelation time of 
several observables O in both models. Since we are comparing two different 
algorithms we normalize the autocorrelation times as in [7] : define one sweep 
as To = 6V4/V configurations, i.e., the average number of attempts needed 
to change every plaquette of the lattice as the unit for the integrated auto- 
correlation times T® nt . In units of updates we have tq = 6V4/ '(V N 'updates)-, 
where N up d ates is defined as either V4 worms for the SWA or 6V4 plaquette 
updates plus 4V4 cube updates for the LMA, i.e., a total of IOV4 updates. 

In order to obtain a measure for the computational effort, the results 
are multiplied by the cost ratio C. In other words we show Tj nt = Crj nt /ro, 
where Tj n j simply is the unnormalized autocorrelation time in units of up- 
dates. The statistical errors of autocorrelation times were estimated with a 
jackknife analysis and were found at the 10 percent level for the statistics 
at our disposal. This is sufficient for the subsequent comparison of the two 
algorithms. 

For the autocorrelation analysis we use the same sets and volumes as 
for the discussion of the characteristic quantities of the SWA and the LMA 
in the previous subsection. Table [3] shows the autocorrelation times in the 
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Figure 10: Normalized histograms of the worm length for the Z3 model 
(upper plot) and the U(l) model (lower plot). 

Z3 case for the SWA and Table [4] is for the LMA. Similarly, Tables [5] and [6] 
correspond to the U(l) case. 

First, we observe that the autocorrelation times for the set close to the 
first order transition (set U-2) increase with the volume, while the others 
are essentially volume independent. It is also interesting to look at the 
sets Z-2 and U-3, where Vlma approaches 6 (see Tables 1 and 2), i.e., the 
configuration space is dominated by closed surfaces, since boundary flux is 
costly for these parameter sets. On the one hand, r^ t and rf^ t are larger 
for the worm algorithm, which is due to the fact that the worm updates 
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Parameters 


V 


int 


int 


int 


int 


oet: Zi-1 


4 X OU 


Lib 


y i 


n n 

u.y 


U.D 


fi, — U.l, 


o x ou 




on 


J..U 


U.D 


ft — n 6 


19 3 x 50 


900 




1 


fi 


/x = 2.0 


16 3 x 50 


197 


88 


1.3 


0.8 


Set: Z-jg 


4 3 x 50 


81 


36 


25 


13 


fc = 0.1, 


8 3 x 50 


84 


32 


25 


14 


P = 0.8, 


12 3 x 50 


> 83 


38 


27 


12 


H = 1.6 


16 3 x 50 


> 90 


37 


30 


13 


Set: Z-3 


4 4 


2.5 


1.3 


0.3 


0.2 


k = 0.5, 


8 4 


5.4 


2.9 


0.6 


0.4 


/3 = 0.28, 


12 4 


6.0 


3.1 


0.6 


0.5 


n = 0.0 


16 4 


7.7 


3.2 


0.7 


0.5 



Table 3: Z3 model: SWA autocorrelation times for different parameter sets. 



Parameters V r u int r? f r*" t 



Set: Z-l 


4 3 x 50 


5429 


2858 


5572 


3096 


K = 0.1, 


8 3 x 50 


5810 


3000 


5858 


3238 


/3 = 0.6, 


12 3 x 50 


5430 


3032 


6072 


4175 


/i = 2.0 


16 3 x 50 


5430 


3048 


>7821 


4256 


Set: Z-2 


4 3 x 50 


67 


48 


748 


306 


K = 0.1, 


8 3 x 50 


68 


51 


753 


302 


= 0.8, 


12 3 x 50 


70 


49 


604 


348 


H = 1.6 


16 3 x 50 


71 


46 


602 


343 


Set: Z-3 


4 4 


107 


55 


59 


23 


k = 0.5, 


8 4 


114 


66 


65 


24 


/3 = 0.28, 


12 4 


116 


69 


67 


25 


n = 0.0 


16 4 


129 


73 


67 


27 



Table 4: Z3 model: LMA autocorrelation times for different parameter sets. 



links in every move, so if the Boltzmann weight of the link variables is very 
low then most of the worms have only a few segments (see Fig. 10). On the 
other hand Ti n t of the observables that depend only on the link occupation 
number is much smaller for the SWA, a fact which reflects the very low 
acceptance rate of the plaquette update of the LMA. 
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Parameters 


V 


={/ 

int 


int 




int 


Set: U-1 


/I 4 
4 


l.Z 


O.l 


U.D 


n q 
U.o 


K = 5, 


«4 
o 


Z.O 


1.0 


U.O 


U.o 


\ — 1 

A — 1, 


1 9 4 

_l_ Z, 






u.u 


u.o 


/3 = 0.40 


16 4 


2.6 


1.1 


0.5 


0.4 


Set: C/"-jS 


4 4 


5.7 


3.5 


9.5 


3.9 


k = 5, 


8 4 


12 


6.9 


2.9 


1.2 


A = l, 


12 4 


19 


7.8 


3.3 


1.4 


/3 = 0.65 


16 4 


21 


7.9 


3.1 


1.6 


Set: C/-5 


4 4 


1584 


865 


1.1 


0.9 


K = 8, 


8 4 


1672 


839 


1.2 


1.0 


A = l, 


12 4 


>1595 


742 


1.8 


0.9 


p = 1.10 


16 4 


>1650 


820 


2.2 


1.1 



Table 5: U(l) model: SWA autocorrelation times for different parameters. 



Parameters 




T U 

int 


int 


int 


T *I0I 2 

int 


Set: U-1 


4 4 


5776 


2926 


8102 


4544 


K = 5, 


8 4 


5809 


3021 


8569 


4532 


A = 1.0, 


12 4 


6122 


4107 


7363 


4959 


/3 = 0.4 


16 4 


6198 


4184 


9065 


5036 


Set: U-2 


4 4 


71 


48 


183 


93 


K = 5, 


8 4 


4737 


2640 


7126 


4108 


A = 1.0, 


12 4 


7167 


2766 


8718 


4275 


/3 = 0.65 


16 4 


7251 


2808 


9430 


4987 


Set: C/-5 


4 4 


464 


282 


446 


301 


K = 8, 


8 4 


428 


302 


480 


294 


A = l, 


12 4 


693 


285 


445 


267 


p = 1.10 


16 4 


711 


280 


489 


271 



Table 6: U(l) model: LMA autocorrelation times for different parameters. 

In general, comparing the results of both algorithms for the two different 
models, we can conclude that the SWA outperforms the LMA for a large 
range of parameters. Only in the region of the space of couplings where the 
link weight is very large the worm algorithm has difficulties to sample the 
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system efficiently, a problem which can easily be overcome with extra cube 
sweeps or by adding a worm with only plaquettes suggested in [JJ. 

5 Summary 

In this article we present a generalization of the worm algorithm to systems 
that are described by surfaces with boundaries of flux, i.e., abelian Gauge- 
Higgs systems. Rewriting the standard form of abelian Gauge-Higgs systems 
in terms of surfaces and fluxes (dual representation) overcomes the complex 
action problem at finite chemical potential. We study Gauge-Higgs systems 
with two gauge groups Z3 and U(l). For the Z3 case a chemical potential 
can be coupled and the system has a complex action problem. 

The key idea of our newly developed surface worm algorithm (SWA) 
is to build up filament-like structures where the dual degrees of freedom 
are changed by adding segments built from plaquette variables and two 
lines of matter flux. We compare the SWA to a local Metropolis algorithm 
(LMA) for the dual representation and in the cases without a sign problem 
also to a conventional Monte Carlo simulation in the standard approach. 
The comparison is used to establish the correctness of the SWA in several 
simulations at different parameter values. 

To study the performance of the SWA we analyze characteristic quan- 
tities: the starting probability, the number of updated plaquettes and the 
cost efficiency. Based on these characteristic quantities we conclude that 
for both gauge groups and most parameter values the SWA is considerably 
more efficient than the LMA. This finding is confirmed by an analysis of 
autocorrelation times where again the SWA is found to decorrelate faster 
(partly considerably faster) than the LMA. 

We expect that the generalization of the worm concept to surface-type 
degrees of freedom will contribute to developing new tools for systems with 
gauge interactions in a dual language. Another important aspect is that 
models where the complex action problem is solved may serve as reference 
systems for testing other approaches such as various reweighting and expan- 
sion techniques. 
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